4. Host reproduction

In this section we explored how the exposure and the infection by heat-stressed parasites affected the host reproduction. First we evaluated the number of individuals fully castrated in exposed uninfected and in infected Daphnias. This divisions into two sub analysis helped us to find more adapted models to describe the host reproductive patern.

4.1 Full castration

One of the consequences of Pasteuria ramosa infection in Daphnia is the reduction in allocation to reproduction, which can also potentially be due to the energetic cost of immune response for uninfected hosts. Here we focused on the proportion of non-reproducing individuals.

4.1.1 Infected individuals

FC1 = glmmTMB(data=infected, Never_repro ~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=infected, Never_repro ~ Isolate + Temperature + DeathAge + (1|Batch), family = "binomial", REML=T)

AICc(FC1, FC1R)
##      df     AICc
## FC1   7 239.8814
## FC1R  8 238.7983
FC1 = glmmTMB(data=infected, Never_repro ~ Isolate * Temperature + DeathAge, family = "binomial", 
               na.action = na.fail)
model_selection = dredge(FC1)

# see best models
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge, 
##     data = infected, family = "binomial", na.action = na.fail, 
##     ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df   logLik
## 16    0.79900          + -0.02228        +        +            + 13  -98.284
## 15   -0.19180          +                 +        +            + 12  -99.622
## 4     0.98420          + -0.01982        +                        5 -108.183
## 3     0.11780          +                 +                        4 -109.410
## 2     0.97250          + -0.02725                                 2 -111.825
## 8     0.94670          + -0.01938        +        +               7 -107.700
## 7     0.09487          +                 +        +               6 -108.859
## 1    -0.25280          +                                          1 -114.432
## 6     0.91880          + -0.02678                 +               4 -111.381
## 5    -0.28880          +                          +               3 -113.875
##     AICc delta weight
## 16 224.9  0.00  0.321
## 15 225.3  0.32  0.273
## 4  226.7  1.79  0.131
## 3  227.1  2.12  0.111
## 2  227.7  2.78  0.080
## 8  230.1  5.16  0.024
## 7  230.2  5.30  0.023
## 1  230.9  5.94  0.016
## 6  231.0  6.06  0.016
## 5  233.9  8.95  0.004
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  Never_repro Never_repro Never_repro
Predictors Odds Ratios CI p Odds Ratios CI p Odds Ratios CI p
DeathAge 0.98 0.95 – 1.00 0.104 0.98 0.96 – 1.01 0.120
Isolate [C1] 2.27 0.66 – 7.83 0.193 2.46 0.72 – 8.38 0.151 1.23 0.49 – 3.10 0.658
Isolate [C18] 0.55 0.18 – 1.69 0.298 0.50 0.16 – 1.49 0.212 0.44 0.18 – 1.11 0.083
Isolate [C24] 0.65 0.20 – 2.09 0.472 0.66 0.21 – 2.08 0.474 0.50 0.20 – 1.28 0.150
Temperature [linear] 0.17 0.03 – 0.96 0.044 0.14 0.03 – 0.77 0.024
Temperature [quadratic] 0.28 0.07 – 1.12 0.073 0.29 0.07 – 1.16 0.080
Isolate [C1] ×
Temperature [linear] 		32.02 3.12 – 328.91 0.004 36.19 3.56 – 367.88 0.002
Isolate [C18] ×
Temperature [linear] 		2.04 0.26 – 16.14 0.497 2.80 0.37 – 21.15 0.319
Isolate [C24] ×
Temperature [linear] 		5.08 0.56 – 45.81 0.147 6.46 0.74 – 56.47 0.092
Isolate [C1] ×
Temperature [quadratic] 		10.98 1.60 – 75.39 0.015 10.14 1.50 – 68.66 0.018
Isolate [C18] ×
Temperature [quadratic] 		3.86 0.64 – 23.40 0.142 3.51 0.59 – 20.92 0.168
Isolate [C24] ×
Temperature [quadratic] 		3.68 0.59 – 23.04 0.163 3.31 0.54 – 20.28 0.196
Observations 167 167 167
R2 Tjur 0.178 0.163 0.074 
best_model = get.models(model_selection, 2)[[1]]
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3977238 0.8536237 0.8579498 0.6620259 0.61213 0.9357635 0.05493602 0.2030844 0.1008657 0.9558345 0.1985479 0.8851658 0.4544661 0.1936561 0.1833452 0.2271874 0.5898748 0.8261575 0.5827213 0.7709696 ...
tab_model(best_model, show.intercept = F)
  Never_repro
Predictors Odds Ratios CI p
Isolate [C1] 2.46 0.72 – 8.38 0.151
Isolate [C18] 0.50 0.16 – 1.49 0.212
Isolate [C24] 0.66 0.21 – 2.08 0.474
Temperature [linear] 0.14 0.03 – 0.77 0.024
Temperature [quadratic] 0.29 0.07 – 1.16 0.080
Isolate [C1] ×
Temperature [linear] 		36.19 3.56 – 367.88 0.002
Isolate [C18] ×
Temperature [linear] 		2.80 0.37 – 21.15 0.319
Isolate [C24] ×
Temperature [linear] 		6.46 0.74 – 56.47 0.092
Isolate [C1] ×
Temperature [quadratic] 		10.14 1.50 – 68.66 0.018
Isolate [C18] ×
Temperature [quadratic] 		3.51 0.59 – 20.92 0.168
Isolate [C24] ×
Temperature [quadratic] 		3.31 0.54 – 20.28 0.196
Observations 167
R2 Tjur 0.163
Visualisation
## `summarise()` has grouped output by 'Temperature'. You can override using the
## `.groups` argument.


A heatwave experienced by the parasite appears to relax its castrating effect on the host for most genotypes, with a strong opposit effect for C1 genotyp, with a lower probability of total sterility in genotype C1, which shows a dramatically increased risk of host sterility after exposure to 40°C. This suggests an interaction between genotype and heatwave experienced by the parasite.

4.1.2 Uninfected individuals

options(na.action = "na.fail")

FC1 = glmmTMB(data=not_infected, Never_repro~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=not_infected, Never_repro ~ Isolate + Temperature+ DeathAge + (1|Batch), family = "binomial", REML=T)

AICc(FC1, FC1R)
##      df     AICc
## FC1   7 118.8705
## FC1R  8 106.9146
FC1 = glmmTMB(data=not_infected, Never_repro ~ Isolate * Temperature + DeathAge + (1|Batch), family = "binomial")
model_selection = dredge(FC1)

# Table of model comparison
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge + 
##     (1 | Batch), data = not_infected, family = "binomial", ziformula = ~0, 
##     dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df  logLik
## 2       1.672          +  -0.1281                                 3 -48.375
## 6       1.459          +  -0.1208                 +               5 -46.911
## 4       2.324          +  -0.1529        +                        6 -46.318
## 8       2.030          +  -0.1347        +        +               8 -45.341
## 16     13.600          +  -0.8717        +        +            + 14 -40.057
## 5      -2.909          +                          +               4 -67.021
## 7      -2.961          +                 +        +               7 -65.982
## 1      -2.830          +                                          2 -72.618
## 15     -3.165          +                 +        +            + 13 -61.318
## 3      -2.927          +                 +                        5 -71.734
##     AICc delta weight
## 2  102.9  0.00  0.499
## 6  104.1  1.24  0.268
## 4  105.0  2.17  0.168
## 8  107.4  4.51  0.052
## 16 110.2  7.32  0.013
## 5  142.2 39.37  0.000
## 7  146.5 43.63  0.000
## 1  149.3 46.43  0.000
## 15 150.4 47.56  0.000
## 3  153.8 50.89  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | Batch)
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3977238 0.8536237 0.8579498 0.6620259 0.61213 0.9357635 0.05493602 0.2030844 0.1008657 0.9558345 0.1985479 0.8851658 0.4544661 0.1936561 0.1833452 0.2271874 0.5898748 0.8261575 0.5827213 0.7709696 ...
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  Never_repro Never_repro
Predictors Odds Ratios CI p Odds Ratios CI p
DeathAge 0.88 0.82 – 0.94 <0.001 0.89 0.83 – 0.95 0.001
Temperature [linear] 0.36 0.10 – 1.30 0.120
Temperature [quadratic] 0.61 0.17 – 2.17 0.442
Random Effects
σ2 3.29 3.29
τ00 6.14 Batch 5.52 Batch
ICC 0.65 0.63
N 50 Batch 50 Batch
Observations 218 218
Marginal R2 / Conditional R2 0.385 / 0.786 0.421 / 0.784 
# best model
best_model = get.models(model_selection, 1)[[1]]
summary(best_model)
##  Family: binomial  ( logit )
## Formula:          Never_repro ~ DeathAge + (1 | Batch)
## Data: not_infected
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.8    112.9    -48.4     96.8      215 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Batch  (Intercept) 6.137    2.477   
## Number of obs: 218, groups:  Batch, 50
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.67209    0.98865   1.691  0.09078 .  
## DeathAge    -0.12808    0.03582  -3.575  0.00035 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(best_model)
  Never_repro
Predictors Odds Ratios CI p
(Intercept) 5.32 0.77 – 36.96 0.091
DeathAge 0.88 0.82 – 0.94 <0.001
Random Effects
σ2 3.29
τ00 Batch 6.14
ICC 0.65
N Batch 50
Observations 218
Marginal R2 / Conditional R2 0.385 / 0.786


Only the age at death seems to reprduce at least once in their life time.

Visualisation
## `summarise()` has grouped output by 'DeathAgeGroup'. You can override using the
## `.groups` argument.

4.1.3 Uninfected individuals (including unexposed husbandry control)

options(na.action = "na.fail")

dataC_ninf = subset(dataC, InfectedStatus!="infected")
summary(dataC_ninf)
##        ID       Temperature Isolate       Clutch        Nb_offspring   
##  Min.   :   2   20: 84      C1  :45   Min.   : 1.000   Min.   :  0.00  
##  1st Qu.: 218   30: 73      C14 :58   1st Qu.: 3.750   1st Qu.: 10.00  
##  Median : 668   40:100      C18 :67   Median : 8.000   Median : 47.00  
##  Mean   : 588               C24 :60   Mean   : 7.368   Mean   : 43.02  
##  3rd Qu.: 969               CTRL:27   3rd Qu.:11.000   3rd Qu.: 68.00  
##  Max.   :1100                         Max.   :16.000   Max.   :111.00  
##                                       NA's   :29                       
##  FirstClutchAge  LastClutchAge      DeathAge     Accidental_death
##  Min.   : 1.00   Min.   :10.00   Min.   :15.00   Min.   :1       
##  1st Qu.:14.00   1st Qu.:39.00   1st Qu.:31.00   1st Qu.:1       
##  Median :18.00   Median :59.00   Median :62.00   Median :1       
##  Mean   :21.77   Mean   :50.47   Mean   :51.47   Mean   :1       
##  3rd Qu.:27.00   3rd Qu.:67.00   3rd Qu.:67.00   3rd Qu.:1       
##  Max.   :67.00   Max.   :67.00   Max.   :67.00   Max.   :1       
##  NA's   :30      NA's   :32                      NA's   :250     
##  Volume_counted InfectedStatus     Number of mature spores
##  Min.   :150    Length:257         Length:257             
##  1st Qu.:150    Class :character   Class :character       
##  Median :150    Mode  :character   Mode  :character       
##  Mean   :152                                              
##  3rd Qu.:150                                              
##  Max.   :200                                              
##  NA's   :12                                               
##  Number of imature spores Nb of white spores cauliflowers  grapes   
##  Length:257               Length:257         no  :236     no  :236  
##  Class :character         Class :character   NA's: 21     NA's: 21  
##  Mode  :character         Mode  :character                          
##                                                                     
##                                                                     
##                                                                     
##                                                                     
##      whites      group      statusDeath     Number of immature spores
##  Min.   :0   C1440  : 29   Min.   :0.0000   Length:257               
##  1st Qu.:0   C2440  : 28   1st Qu.:0.0000   Class :character         
##  Median :0   NACTRL : 27   Median :1.0000   Mode  :character         
##  Mean   :0   C1840  : 26   Mean   :0.5136                            
##  3rd Qu.:0   C1830  : 24   3rd Qu.:1.0000                            
##  Max.   :0   C130   : 19   Max.   :1.0000                            
##              (Other):104                                             
##   Never_repro         Batch          imm         mat     mature_load
##  Min.   :0.0000   2      :  8   Min.   :0   Min.   :0   Min.   :0   
##  1st Qu.:0.0000   17     :  8   1st Qu.:0   1st Qu.:0   1st Qu.:0   
##  Median :0.0000   28     :  8   Median :0   Median :0   Median :0   
##  Mean   :0.1128   31     :  8   Mean   :0   Mean   :0   Mean   :0   
##  3rd Qu.:0.0000   3      :  7   3rd Qu.:0   3rd Qu.:0   3rd Qu.:0   
##  Max.   :1.0000   6      :  7   Max.   :0   Max.   :0   Max.   :0   
##                   (Other):211                           NA's   :12  
##     imm_load      ratio       repro_rate     stade_max 
##  Min.   :0    Min.   : NA   Min.   :0.0000   None:257  
##  1st Qu.:0    1st Qu.: NA   1st Qu.:0.3095             
##  Median :0    Median : NA   Median :0.8358             
##  Mean   :0    Mean   :NaN   Mean   :0.7206             
##  3rd Qu.:0    3rd Qu.: NA   3rd Qu.:1.0597             
##  Max.   :0    Max.   : NA   Max.   :1.6567             
##  NA's   :12   NA's   :257
FC1 = glmmTMB(data=dataC_ninf, Never_repro~ Isolate + Temperature + DeathAge, family = "binomial", REML=T)
FC1R = glmmTMB(data=dataC_ninf, Never_repro ~ Isolate + Temperature+ DeathAge + (1|Batch), family = "binomial", REML=T)

AICc(FC1, FC1R)
##      df     AICc
## FC1   8 134.3687
## FC1R  9 125.3196
FC1 = glmmTMB(data=not_infected, Never_repro ~ Isolate * Temperature + DeathAge + (1|Batch), family = "binomial")
model_selection = dredge(FC1)

# Table of model comparison
model_selection
## Global model call: glmmTMB(formula = Never_repro ~ Isolate * Temperature + DeathAge + 
##     (1 | Batch), data = not_infected, family = "binomial", ziformula = ~0, 
##     dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df  logLik
## 2       1.672          +  -0.1281                                 3 -48.375
## 6       1.459          +  -0.1208                 +               5 -46.911
## 4       2.324          +  -0.1529        +                        6 -46.318
## 8       2.030          +  -0.1347        +        +               8 -45.341
## 16     13.600          +  -0.8717        +        +            + 14 -40.057
## 5      -2.909          +                          +               4 -67.021
## 7      -2.961          +                 +        +               7 -65.982
## 1      -2.830          +                                          2 -72.618
## 15     -3.165          +                 +        +            + 13 -61.318
## 3      -2.927          +                 +                        5 -71.734
##     AICc delta weight
## 2  102.9  0.00  0.499
## 6  104.1  1.24  0.268
## 4  105.0  2.17  0.168
## 8  107.4  4.51  0.052
## 16 110.2  7.32  0.013
## 5  142.2 39.37  0.000
## 7  146.5 43.63  0.000
## 1  149.3 46.43  0.000
## 15 150.4 47.56  0.000
## 3  153.8 50.89  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | Batch)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  Never_repro Never_repro
Predictors Odds Ratios CI p Odds Ratios CI p
DeathAge 0.88 0.82 – 0.94 <0.001 0.89 0.83 – 0.95 0.001
Temperature [linear] 0.36 0.10 – 1.30 0.120
Temperature [quadratic] 0.61 0.17 – 2.17 0.442
Random Effects
σ2 3.29 3.29
τ00 6.14 Batch 5.52 Batch
ICC 0.65 0.63
N 50 Batch 50 Batch
Observations 218 218
Marginal R2 / Conditional R2 0.385 / 0.786 0.421 / 0.784 
best_model = get.models(model_selection, 1)[[1]]
summary(best_model)
##  Family: binomial  ( logit )
## Formula:          Never_repro ~ DeathAge + (1 | Batch)
## Data: not_infected
## 
##      AIC      BIC   logLik deviance df.resid 
##    102.8    112.9    -48.4     96.8      215 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  Batch  (Intercept) 6.137    2.477   
## Number of obs: 218, groups:  Batch, 50
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.67209    0.98865   1.691  0.09078 .  
## DeathAge    -0.12808    0.03582  -3.575  0.00035 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
tab_model(best_model)
  Never_repro
Predictors Odds Ratios CI p
(Intercept) 5.32 0.77 – 36.96 0.091
DeathAge 0.88 0.82 – 0.94 <0.001
Random Effects
σ2 3.29
τ00 Batch 6.14
ICC 0.65
N Batch 50
Observations 218
Marginal R2 / Conditional R2 0.385 / 0.786 
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.9088821 0.9008778 0.9317442 0.7483929 0.241452 0.2361889 0.7475164 0.4387352 0.390151 0.5425091 0.609537 0.7492443 0.6104253 0.7648734 0.5350931 0.2879339 0.134203 0.1612245 0.3157523 0.1570396 ...


Only the age at death seems to reprduce at least once in their life time, the husbandry control inclusion does not change this.

Visualisation
## `summarise()` has grouped output by 'DeathAgeGroup'. You can override using the
## `.groups` argument.

4.2 Number of individuals produced

In this section we were interested in understanding the role of parasite heat-stress and genotype on the number of offspring produced. Indeed, one of the consequences of infection in Daphnia is the reduction in allocation to parasite reproduction, (or a potential cost of immune response for uninfected individuals), even if the host is not fully castrated.

4.2.1 Infected individuals

repro_inf = subset(data, data$Clutch>0 & InfectedStatus=="infected" & Temperature!="NA")
repro_inf$Temperature = factor(repro_inf$Temperature, ordered = T)

FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "nbinom2", REML=T)
#FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "poisson", REML=T)
#FC1 = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "gaussian", REML=T)

simulateResiduals(FC1, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.8741857 0.7996104 0.4441007 0.4213267 0.3856629 0.2460142 0.4371622 0.9950848 0.5756341 0.3335869 0.1840254 0.6387504 0.2376711 0.8736609 0.9228149 0.7387757 0.2361033 0.4922073 0.4105355 0.1956212 ...
FC1R = glmmTMB(data=repro_inf, Nb_offspring ~ Isolate + Temperature + DeathAge + (1|Batch), family = "nbinom2", REML=T)
AICc(FC1, FC1R)# no batch effect
##      df     AICc
## FC1   8 601.8980
## FC1R  9 601.5721
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge  , 
                data = repro_inf, family = "nbinom2", na.action = na.fail)
model_selection = dredge(FC1)

model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge, 
##     data = repro_inf, family = "nbinom2", na.action = na.fail, 
##     ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df   logLik
## 4      0.9519          +  0.03153        +                        6 -283.958
## 2      0.2723          +  0.03680                                 3 -289.038
## 8      0.9455          +  0.03176        +        +               8 -283.902
## 16     0.7038          +  0.02896        +        +            + 14 -277.332
## 6      0.3050          +  0.03628                 +               5 -288.741
## 3      2.7120          +                 +                        5 -291.710
## 7      2.7390          +                 +        +               7 -291.682
## 15     2.2290          +                 +        +            + 13 -284.113
## 1      2.1350          +                                          2 -300.042
## 5      2.1400          +                          +               4 -299.296
##     AICc delta weight
## 4  580.9  0.00  0.750
## 2  584.3  3.46  0.133
## 8  585.5  4.62  0.075
## 16 588.0  7.10  0.022
## 6  588.2  7.28  0.020
## 3  594.1 13.22  0.001
## 7  598.7 17.79  0.000
## 15 598.8 17.89  0.000
## 1  604.2 23.34  0.000
## 5  607.0 26.16  0.000
## Models ranked by AICc(x)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
tab_model(selected_models, show.intercept = F)
  Nb_offspring
Predictors Incidence Rate Ratios CI p
DeathAge 1.03 1.02 – 1.05 <0.001
Isolate [C1] 0.44 0.22 – 0.89 0.021
Isolate [C18] 0.48 0.28 – 0.82 0.008
Isolate [C24] 0.80 0.45 – 1.39 0.424
Observations 94 
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T) 

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.8445634 0.7757172 0.3188654 0.4678265 0.4046863 0.2574613 0.3273075 0.988557 0.6029993 0.2770807 0.2222652 0.5937299 0.200338 0.8887268 0.8925955 0.6897923 0.2627872 0.457771 0.3542096 0.1563993 ...
tab_model(best_model, show.intercept = F)
  Nb_offspring
Predictors Incidence Rate Ratios CI p
DeathAge 1.03 1.02 – 1.05 <0.001
Isolate [C1] 0.44 0.22 – 0.89 0.021
Isolate [C18] 0.48 0.28 – 0.82 0.008
Isolate [C24] 0.80 0.45 – 1.39 0.424
Observations 94 
summary(best_model)
##  Family: nbinom2  ( log )
## Formula:          Nb_offspring ~ DeathAge + Isolate + 1
## Data: repro_inf
## 
##      AIC      BIC   logLik deviance df.resid 
##    579.9    595.2   -284.0    567.9       88 
## 
## 
## Dispersion parameter for nbinom2 family ():  1.4 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.951947   0.468983   2.030  0.04238 *  
## DeathAge     0.031525   0.007622   4.136 3.54e-05 ***
## IsolateC1   -0.812114   0.352549  -2.304  0.02125 *  
## IsolateC18  -0.738421   0.276490  -2.671  0.00757 ** 
## IsolateC24  -0.228773   0.286424  -0.799  0.42445    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Visualisation

4.2.2 Uninfected individuals

repro_nf = subset(data, data$Clutch>0 & InfectedStatus=="not_inf" & Temperature!="NA")
repro_nf$Temperature = factor(repro_nf$Temperature, ordered = T)

FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature + DeathAge+ (1|Batch), family = "gaussian", REML=T)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature+ DeathAge + (1|Batch), family = "poisson")
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate * Temperature + DeathAge+ (1|Batch), family = "nbinom2", REML=T)
#FC1 = glmmTMB(data=repro_nf, log(Nb_offspring) ~ Isolate * Temperature + DeathAge + (1|Batch), family = "gaussian", REML=T)

FC1R = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature + DeathAge  + (1|Batch), family = "gaussian", REML=F)
FC1 = glmmTMB(Nb_offspring ~ Isolate + Temperature + DeathAge, data = repro_nf, family = "gaussian", REML=F)

AICc(FC1, FC1R)# Batch effect NOT important! But significantly improved the residuals
##      df     AICc
## FC1   8 1659.382
## FC1R  9 1660.230
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge + (1|Batch), data = repro_nf, family = "gaussian", na.action = na.fail)
model_selection = dredge(FC1)

model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge + 
##     (1 | Batch), data = repro_nf, family = "gaussian", na.action = na.fail, 
##     ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df   logLik
## 2      -22.63          +    1.306                                 4 -823.741
## 4      -20.37          +    1.313        +                        7 -821.201
## 6      -22.45          +    1.294                 +               6 -823.160
## 8      -20.41          +    1.303        +        +               9 -820.623
## 5       48.57          +                          +               5 -913.919
## 1       50.26          +                                          3 -916.804
## 7       51.89          +                 +        +               8 -913.084
## 3       53.77          +                 +                        6 -915.627
## 15      52.47          +                 +        +            + 14 -912.123
## 16     -20.73          +    1.307        +        +            + 15         
##      AICc  delta
## 2  1655.7   0.00
## 4  1657.0   1.31
## 6  1658.8   3.08
## 8  1660.2   4.54
## 5  1838.2 182.46
## 1  1839.7 184.04
## 7  1843.0 187.26
## 3  1843.7 188.01
## 15 1854.6 198.91
## 16              
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | Batch)
models_within_2_AIC = subset(model_selection, delta < 2)
selected_models = get.models(models_within_2_AIC, 1:nrow(models_within_2_AIC))
selected_models
## $`2`
## Formula:          Nb_offspring ~ DeathAge + (1 | Batch)
## Data: repro_nf
##       AIC       BIC    logLik  df.resid 
## 1655.4818 1668.5326 -823.7409       189 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups   Name        Std.Dev.
##  Batch    (Intercept)  4.599  
##  Residual             16.721  
## 
## Number of obs: 193 / Conditional model: Batch, 50
## 
## Dispersion estimate for gaussian family (sigma^2):  280 
## 
## Fixed Effects:
## 
## Conditional model:
## (Intercept)     DeathAge  
##     -22.631        1.306  
## 
## $`4`
## Formula:          Nb_offspring ~ DeathAge + Isolate + (1 | Batch)
## Data: repro_nf
##       AIC       BIC    logLik  df.resid 
## 1656.4026 1679.2414 -821.2013       186 
## Random-effects (co)variances:
## 
## Conditional model:
##  Groups   Name        Std.Dev.
##  Batch    (Intercept)  4.633  
##  Residual             16.482  
## 
## Number of obs: 193 / Conditional model: Batch, 50
## 
## Dispersion estimate for gaussian family (sigma^2):  272 
## 
## Fixed Effects:
## 
## Conditional model:
## (Intercept)     DeathAge    IsolateC1   IsolateC18   IsolateC24  
##    -20.3723       1.3133      -4.2733      -0.6133      -6.5511  
## 
## $<NA>
## NULL
## 
## attr(,"rank")
## function (x) 
## do.call("rank", list(x))
## <environment: 0x137323810>
## attr(,"call")
## AICc(x)
## attr(,"class")
## [1] "function"     "rankFunction"
## attr(,"beta")
## [1] "none"
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.78 0.656 0.584 0.568 0.468 0.564 0.596 0.448 0.512 0.388 0.464 0.368 0.504 0.452 0.672 0.536 0.708 0.592 0.356 0.488 ...
tab_model(best_model, show.intercept = F) # No effect
  Nb_offspring
Predictors Estimates CI p
DeathAge 1.31 1.16 – 1.45 <0.001
Random Effects
σ2 279.59
τ00 Batch 21.15
ICC 0.07
N Batch 50
Observations 193
Marginal R2 / Conditional R2 0.622 / 0.649 
summary(best_model)
##  Family: gaussian  ( identity )
## Formula:          Nb_offspring ~ DeathAge + (1 | Batch)
## Data: repro_nf
## 
##      AIC      BIC   logLik deviance df.resid 
##   1655.5   1668.5   -823.7   1647.5      189 
## 
## Random effects:
## 
## Conditional model:
##  Groups   Name        Variance Std.Dev.
##  Batch    (Intercept)  21.15    4.599  
##  Residual             279.59   16.721  
## Number of obs: 193, groups:  Batch, 50
## 
## Dispersion estimate for gaussian family (sigma^2):  280 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -22.63107    4.34365   -5.21 1.89e-07 ***
## DeathAge      1.30586    0.07379   17.70  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Calculation of daily reproduction rate
repro_nf$daily_repro_rate = repro_nf$Nb_offspring / repro_nf$DeathAge
mean(repro_nf$daily_repro_rate, na.rm = TRUE)
## [1] 0.8219201
(sem = sd(repro_nf$daily_repro_rate, na.rm = TRUE) / sqrt(length(repro_nf$daily_repro_rate)))
## [1] 0.02774078
Visualisation

Only the effect of time is important here; each day, a Daphnia produces on average 0.30 more offspring.

4.2.3 Uninfected individuals (including unexposed husbandry control)

dataC$Isolate = as.character(dataC$Isolate)
dataC$Isolate[is.na(dataC$Isolate)] = "unexposed"
dataC$Isolate = as.factor(dataC$Isolate)

dataC$Temperature = as.character(dataC$Temperature)
dataC$Temperature[dataC$Temperature == "CTRL"] = 20
dataC$Temperature = factor(dataC$Temperature, ordered = TRUE)

repro_nfC = subset(dataC, dataC$Clutch>0 & InfectedStatus!="inf")

repro_nfC$Temperature = factor(repro_nfC$Temperature, labels = c("20", "30","40"), ordered = T)

FC1 = glmmTMB(data=repro_nfC, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "gaussian", REML=F)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature+ DeathAge, family = "poisson", REML=F)
#FC1 = glmmTMB(data=repro_nf, Nb_offspring ~ Isolate + Temperature + DeathAge, family = "nbinom2", REML=T)

simulateResiduals(FC1, plot=T) # I don't find better and I have no reason to exclude the outlier.

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.724 0.608 0.58 0.488 0.6 0.48 0.656 0.54 0.672 0.492 0.484 0.408 0.4 0.548 0.464 0.352 0.8 0.624 0.824 0.668 ...
FC1R = glmmTMB(data=repro_nfC, Nb_offspring ~ Isolate + Temperature + DeathAge  + (1|Batch), family = "gaussian", REML=F)
FC1 = glmmTMB(Nb_offspring ~ Isolate + Temperature + DeathAge, data = repro_nfC, family = "gaussian", REML=F)

AICc(FC1, FC1R)# Batch effect NOT important and do not improve residuals anyway.
##      df     AICc
## FC1   9 1945.748
## FC1R 10 1946.820
FC1 = glmmTMB(Nb_offspring ~ Isolate * Temperature + DeathAge, data = repro_nfC, family = "gaussian", na.action = na.fail)
model_selection = dredge(FC1) 

model_selection
## Global model call: glmmTMB(formula = Nb_offspring ~ Isolate * Temperature + DeathAge, 
##     data = repro_nfC, family = "gaussian", na.action = na.fail, 
##     ziformula = ~0, dispformula = ~1)
## ---
## Model selection table 
##    cnd((Int)) dsp((Int)) cnd(DtA) cnd(Isl) cnd(Tmp) cnd(Isl:Tmp) df    logLik
## 2      -22.74          +    1.305                                 3  -966.679
## 4      -24.62          +    1.311        +                        7  -963.982
## 6      -22.39          +    1.295                 +               5  -966.198
## 8      -24.57          +    1.305        +        +               9  -963.461
## 16     -24.28          +    1.314        +        +            + 15  -962.257
## 5       47.93          +                          +               4 -1078.562
## 1       48.49          +                                          2 -1083.107
## 7       50.74          +                 +        +               8 -1077.345
## 3       52.70          +                 +                        6 -1079.491
## 15      46.80          +                 +        +            + 14 -1076.437
##      AICc  delta weight
## 2  1939.5   0.00  0.682
## 4  1942.5   3.01  0.151
## 6  1942.7   3.20  0.137
## 8  1945.7   6.28  0.029
## 16 1956.8  17.31  0.000
## 5  2165.3 225.84  0.000
## 1  2170.3 230.80  0.000
## 7  2171.3 231.88  0.000
## 3  2171.4 231.90  0.000
## 15 2182.8 243.38  0.000
## Models ranked by AICc(x)
best_model = get.models(model_selection, 1)[[1]]
simulateResiduals(best_model, plot=T)

## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.776 0.628 0.604 0.564 0.528 0.548 0.648 0.632 0.548 0.436 0.528 0.432 0.496 0.436 0.432 0.452 0.724 0.516 0.784 0.664 ...
selected_models = get.models(model_selection, 1:nrow(model_selection))
tab_model(selected_models, show.intercept = F) # No effect
  Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring Nb_offspring
Predictors Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p Estimates CI p
DeathAge 1.30 1.18 – 1.43 <0.001 1.31 1.18 – 1.44 <0.001 1.29 1.17 – 1.42 <0.001 1.31 1.18 – 1.43 <0.001 1.31 1.18 – 1.44 <0.001
Isolate [C14] 4.03 -2.97 – 11.04 0.259 3.83 -3.22 – 10.88 0.287 3.21 -4.95 – 11.37 0.441 1.19 -10.42 – 12.80 0.841 0.82 -10.79 – 12.43 0.890 5.39 -8.07 – 18.85 0.433
Isolate [C18] 3.72 -3.09 – 10.52 0.284 3.82 -3.00 – 10.63 0.272 3.17 -4.64 – 10.98 0.426 -3.36 -14.52 – 7.80 0.555 -4.35 -15.56 – 6.87 0.447 0.98 -11.90 – 13.86 0.882
Isolate [C24] -2.42 -9.45 – 4.62 0.501 -2.51 -9.64 – 4.63 0.491 -2.96 -11.08 – 5.16 0.475 -5.37 -17.12 – 6.37 0.370 -6.47 -18.12 – 5.19 0.277 -1.23 -14.62 – 12.17 0.858
Isolate [CTRL] 1.21 -7.28 – 9.69 0.780 2.25 -7.51 – 12.00 0.652 -0.81 -18.17 – 16.54 0.927 -7.95 -23.93 – 8.03 0.330 -15.55 -29.36 – -1.74 0.027 7.65 -20.95 – 36.26 0.600
Temperature [linear] 1.63 -2.10 – 5.37 0.391 1.71 -2.54 – 5.96 0.431 -1.60 -14.45 – 11.26 0.808 9.22 3.23 – 15.20 0.003 7.35 0.41 – 14.29 0.038 19.74 -1.19 – 40.67 0.065
Temperature [quadratic] 0.81 -3.19 – 4.80 0.692 0.98 -3.23 – 5.18 0.649 1.01 -8.88 – 10.90 0.841 -2.19 -8.71 – 4.33 0.510 -1.09 -8.01 – 5.83 0.757 -8.20 -24.46 – 8.05 0.323
Isolate [C14] ×
Temperature [linear] 		3.67 -11.49 – 18.83 0.635 -14.13 -38.99 – 10.72 0.265
Isolate [C18] ×
Temperature [linear] 		3.79 -11.12 – 18.70 0.618 -14.78 -39.20 – 9.63 0.235
Isolate [C24] ×
Temperature [linear] 		3.84 -10.96 – 18.64 0.611 -12.61 -36.88 – 11.67 0.309
Isolate [C14] ×
Temperature [quadratic] 		-2.63 -15.81 – 10.55 0.696 8.53 -13.14 – 30.20 0.440
Isolate [C18] ×
Temperature [quadratic] 		4.00 -8.12 – 16.11 0.518 9.27 -10.71 – 29.24 0.363
Isolate [C24] ×
Temperature [quadratic] 		-2.94 -16.34 – 10.46 0.667 5.90 -16.16 – 27.96 0.600
Observations 228 228 228 228 228 228 228 228 228 228
R2 / R2 adjusted 0.640 / 0.637 0.648 / 0.639 0.641 / 0.635 0.650 / 0.637 0.654 / 0.631 0.039 / 0.026 0.000 / -0.004 0.049 / 0.019 0.031 / 0.009 0.057 / -0.000 
summary(best_model)
##  Family: gaussian  ( identity )
## Formula:          Nb_offspring ~ DeathAge + 1
## Data: repro_nfC
## 
##      AIC      BIC   logLik deviance df.resid 
##   1939.4   1949.6   -966.7   1933.4      225 
## 
## 
## Dispersion estimate for gaussian family (sigma^2):  282 
## 
## Conditional model:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -22.74119    3.70948  -6.131 8.76e-10 ***
## DeathAge      1.30473    0.06482  20.127  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

  Again, the husbandry control do not change the conclusions.

Visualisation

4.2.4 Descriptive stat per group

Here we just described the average number of offspring produced per day in total in each group

# Function to compute mean and 95% confidence interval
mean_ci = function(x) {
  x = x[!is.na(x)]
  m = mean(x)
  se = sd(x) / sqrt(length(x))
  c(mean = m,
    CI_low = m - 1.96 * se,
    CI_high = m + 1.96 * se)
}

# Define groups
g1 = mean_ci(repro_nfC$repro_rate[repro_nfC$Isolate == "CTRL"])
g2 = mean_ci(repro_nf$repro_rate[repro_nfC$Temperature == 20])
g3 = mean_ci(repro_nf$repro_rate)
g4 = mean_ci(repro_inf$repro_rate)

rbind(
  noninfected_unexposed     = g1,
  noninfected_exposed_20    = g2,
  noninfected_exposed_all   = g3,
  infected_all              = g4
)
##                              mean    CI_low   CI_high
## noninfected_unexposed   0.7203259 0.5747725 0.8658793
## noninfected_exposed_20  0.7463513 0.6490932 0.8436094
## noninfected_exposed_all 0.8219201 0.7675482 0.8762921
## infected_all            0.1699311 0.1302886 0.2095736